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I. INTRODUCTION 



A . Background 

The Northrop Corporation, under Air Force funding, has developed 
a finite element digital computer code, called BR-1, for predicting the 
inelastic, large deflection, transient response of combat aircraft skin- 
rib-stringer structures when subjected to internal air blast loading. 

The finite elements considered are flat rectangular plates and beam 
stiffeners. The theory, user’s manual and code listing are given in 
References 1 and 2. The Air Force Flight Dynamics Laboratory wanted 
the BR-1 code modified so that it could be used to predict the response 
of aircraft fuel tank walls when subjected to fluid pressures due to 
projectiles passing through the fuel in the tank. The intense pressure 
and momentum in the fuel due to the penetrating projectile is referred 
to as the hydraulic ram loading. This report describes the modifica- 
tions to the IBM version of the BR-1 code to account for the fluid 
(fuel) - structure (tank wall) interaction that occurs when bullets and 
metal fragments penetrate into aircraft fuel tanks. The modified code 
is called BR-1HR. The interaction between the compressible fluid and 
the structure is approximated by the piston theory. The code can also 
be used for many other compressible fluid-structure interaction problems. 

B. Piston Theory 

The total nonlinear problem of the response of a tank containing a 
fluid and subjected to a high speed penetrating projectile is extremely 
complex and presently defies analytical treatment. In general, the 
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equations for the fluid stresses and motion are coupled to those for the 
wall stresses and motion due to the continuity at the fluid- structure 
interface (3). One procedure for approximating the fluid-structure in- 
teraction phenomenon is the piston theory (4) . This theory has been in 
use since the early 1940' s when it was applied to the study of the 
effect of underwater explosions on ship plates. It provides the correct 
solution to the one -dimensional propagation of stresses in an acoustic 
medium due to a moving boundary. Several recent studies have been made 
to determine its accuracy when applied to two dimensional fluid-structure 
interaction problems (4,5). 

Application of the piston theory to the interaction problem allows 
the structure equations and fluid equations to be uncoupled. The response 
of the wall is computed using the conventional structural response equa- 
tions, with the normal pressure on the wall p given by 

p = P o + pc (v ± - w) (1) 

where p and v are the incident pressure and velocity of the fluid at 
o i 

the wall respectively, P is the fluid density, c is the acoustic velocity 

in the fluid, and w is the wall velocity*. The pressure, p^, and 

velocity, v^, are the pressure and velocity that would exist in the 

fluid if the interface was not there, i.e., p and v. do not contain 

o i 

any "local" reflected effects. However, effects on p and v. due to 

o 1 

earlier reflections from other walls and free surfaces should be con- 
sidered. In other words, p Q and v^ are the loading components due to 

*A dot above a variable denotes a derivative with respect to time. 
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the free field and the scattered effects. The loading component due to 

the wall velocity w is called the radiation pressure. 

C. The NWC Hydraulic Ram Computer Code 

In order to use the piston theory to compute the tank wall response, 

it is necessary to know the incident fluid pressure p and velocity v. 

o 1 

over the entire fluid-wall interface as a function of time. In conjunc- 
tion with this project Lundstrom, at the Naval Weapons Center, has 
developed a digital computer code that predicts the fluid pressures and 
velocities p Q and v^ throughout a rectangular body of fluid due to a 
penetrating ballistic projectile. The model is based upon replacing 
the projectile by a line of sources whose strength is determined by an 
energy balance between the kinetic and potential energy of the fluid and 
the energy loss due to drag forces on the projectile. Reflections from 
the structure -fluid interface are accounted for by considering the fluid 
boundary to be either stress free or rigid*. An extensive series of 
tests were performed at the Naval Weapons Center to obtain detailed 
pressure measurements for a variety of projectiles under a wide range of 
impact conditions. This data allowed the selection of important para- 
meters such as tumbling distance, jacket stripping, etc., to be entered 
into the code. A description of the code and the instructions for 
operation are given in Reference 7* This code provides the values for 

p and v. at user specified locations over the structure -fluid interface 
o 1 

for the time span of interest. 

* A study of the one-dimensional reflection of step pressure waves from 
typical aircraft fuel tank walls indicates that the stress free surface 
provides the more accurate approximation (6) 
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II. MODIFICATION OF THE BR-1 CODE 



A. Incorporation of the Piston Theory 

The BR-1 code has an option for the user to input a time varying 
pressure on each panel element. In the piston theory this pressure is 

the p Q + pcv^ of Eq. 1. The other contributor to the wall loading given 

« 

by Eq. 1 is w, the wall velocity. Since the BR-1 code does not include 
damping effects, it is necessary to add the damping term pew to the 
equations of motion. 

The BR-1 code solves the set of equations (Eqs. 1 and 2, Ref. l) 

[M] {q*} = {F} - {P} - [H] {q*} = {C} (2) 

for the vector of global nodal generalized displacements {q*} as a function 
of time. These generalized displacements define the motion of the walls. 
The vector [F] consists of global generalized external and body forces 
at the nodes of the elements. The matrix [M] is the mass matrix. 

The wall pressure p given by Eq. 1 causes external forces at the 
nodes. The external generalized forces at the nodes of each element in 
the local coordinate system, {f}, is given by (Eq. A-47, Ref. l) 




T 

where ft] is the transpose of the matrix of shape functions [N] , 
evaluated at the surface of the element, Aout is the surface of the ele- 
ment, and {T q } is the vector of applied surface tractions and moments. 
The order of {T q } is a 5x1 vector. Due to the fluid pressure loading 
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where the subscripts denote the coordinates, (T is normal to the ele- 
ment), and p is given by Eq. 1. The fourth and fifth elements of {T^} 
correspond to applied moments per unit middle surface area of the 
element, and are zero here. The numerical intergration of Eq. 3 can 
be accomplished by Gaussian quadrature. However, the [f] is obtained 
in the BR-1 code in a more approximate way by using a lumping approach 
at the nodes of the element, as is done for the mass matrix evaluation, 
in order to save computation time. Thus, according to Eq. B-92 of 
Ref. 1, the external force vector at the rth node of the nth element 
is given by 



= ( X 2- X l) ( y 2- y l) ,n -n 

l nri 4 



nr J 



D 



nr 



1 

0 

0 



(5) 



nr 



where p^ is the magnitude of the pressure on the element*, and 



"\/( 1+0 ) nr where 0^ and are the fourth and fifth local general- 
ized displacements at the rth node. They appear here because the 
pressure is defined in BR-1 as the pressure normal to the deformed 
surface. The quantities (x^-x^) and (y -y ) are the dimensions of the 
rectangular element. 



The pressure p^ in the piston theory is given by Eq. 1, i.e. 



p = P + (pc) v. - (pc) w 
n on K n in K n nr 



( 6 ) 



* The assumption is made in the programming of BR-1 that the pressure 
is uniform over each element. This is contrary to the theoretical 
presentation where the pressure is defined at each node point. 
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where 



! pc if the nth element is in contact with the fluid 

0 if the nth element is not in contact with the fluid. 
The variables p and v. can be determined from the NWC computer code 
for each element for the time span of interest prior to the computation 
of the wall response. This data is then input as the known external 
pressure. The variable w is an unknown dependent variable and is part 
of {q*}. Hence, it must be incorporated into the equations of motion, 
Eq. 2. 

The f f 1 due to w is given by 
i nr J nr 




The global force vector at the rth node, {F}^, is related to {f } in 
the form 

fF] = E [J ] T ff } (7b) 

where [J ] is the transformation matrix from the global coordinate 
system to the nth element local coordinate system, g means a summation 
over all elements containing the node r, and a — * r means node <y corre- 
sponds physically to node r. Since the .wall velocity in Eq. 7a is given 
in terms of the local coordinates, it must be converted to global 
coordinates. Thus, according to Eqs. A-77 and B-3 of Reference 1, 

w = fq*} (8) 

nr L n J 1 ; r 
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3 

where L^ n J the third row of [J^]. Thus, the global external force 
vector at the rth node {F}^ becomes 




• • • 

Note that {F}_^ is nonlinear since and 0^ are part of If the 

rotations and ©^ are neglected in the computation of {F}_^, i.e., if 
the pressure is not truly normal to the deformed surface, [F] for the 
total R nodes of the structure can be given in the form 



{F} =. - [D] {4*) 



(10a) 



where 



[D] = 









and [D] r is a 6x6 matrix given by 



m 



(10b) 



P>] r - i S ( p c) ( X 2- X l) ( y 2- y l) n [j ] x LJ n °J 



,T , T 3 



(10c) 
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B. Method of Solution 

The BR-1 code solves for {q*} at discrete points in time using the 
explicit finite difference scheme (Eq. A-109, Ref. l) 



{Aq*} t = {Aq*} t + At {q*} t + (At) (q*} t 

i+1 i i i 

where At is the time interval between two time points, i.e. 

At = t. +1 - t x 



and 



{Aq*} = {q*j - {q*) 

i+1 l+l j 



The acceleration {‘q*}^ is obtained from Eq. (2) in the form 

i 



( 11 ) 



( 12 ) 



{q*} t> = CM]' [C} t _ 



(13) 



The {q*} are due to impulsive loads which are known in the blast loading 

problem. In the BR-1 code {F}, {P}, {q*} and {q|:} are known at time t^. 

Hence, {Aq*}, and {q*} can be determined using Eqs. 11-13. For 
i+1 i+1 

our situation, {F}^ contains {q*}^_ , which is unknown. We could 

i i 

approximate {q*}^_ with the backward finite difference form 



{q*} + = {Aq*}, /At 



t . 
1 



(14) 



If we do, then {q*} becomes known at t,, {F}, and hence { C} , can be 

w • 1 

1 

determined at t^, and the procedure used in BR-1 is directly applicable. 

On the other hand, if we express {q*}^ In the central finite difference 

i 

form 



{q*} t = ( {Aq* } , 

i \ 1+1 



{Aq*} 




(2At) 



(15) 
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. Consequently 



then {F} , and hence {0} , depends upon {Aq] 

i i 1+1 

{Aq*} appears on both the left and right hand side of Eq. 11. This 
l+l 

requires a new solution procedure. A detailed study of the accuracy 
and numerical stability of these two approaches, and a third approach, 
when applied to a single degree of freedom, damped oscillator is pre- 
sented in the Appendix. The approach where {q*} is given by the central 
difference expression, Eq. 15, is the one selected based upon the 
accuracy and stability properties of this scheme. Its shown in the 
Appendix that the maximum value of At for a stable solution is 2/^, where 
uj is the highest natural undamped frequency. This is identical to the 
stability limit on the BR-1 procedure. 

Introducting that part of {F} due to w given by Eq. 10a into Eq. 2 
results in the modified equations of motion 

[M] {q*} + [D] {q*3 = {C} (16) 

Replacing {q*} with the conventional central difference approximation 



{•4*3 = [ {q*3 -2 {q*3 + {q*3 ) /(At) 2 (17a) 

i \ l+l z i l+l )/ 

g f * 

is equivalent to obtaining {q*}^ from Eq. 11 with {q*3 not considered, i.e. 



{q*3 t> = ^ " {^3 t . |^/(At) 2 (17b) 



according to Eq. 12. Substituting Eq. 15 for {q*}^ and Eq. 17b for 



{*4*3+ into Eq. 16 leads to 



{Aq*3 t ^ = [M + D (At/2)]" 1 ^ [M - D (At/2)] {Aq*) t> + {C) t J(l8) 
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which is equivalent to fAq*} given by Eqs. 11 and 13 when D is a 

l+l 

null matrix and {q*} is not considered. 

The mass matrix [M] is developed in BR-1 using the lumped mass 
approach and is given by (Eq. A-97, Ref. l) 



[Ml = 




0 




where [M] 



r 



is a 6x6 matrix given by 




r 



(19a) 



(19b) 



and [m ] is a diagonal matrix of the lumped mass at the a node of the 
n a 

nth element. Comparing Eq. 10b with Eq. 19a reveals that the two 
matrices [M + D (At/2)] and [M - D (At/2)] occupy the same nonzero 6x6 
locations as the original matrix M. Thus, the same procedure used in 
BR-1 to compute [M] ^ can be used to compute [M + D (At/2)] \ Its 
only necessary to modify the elements of [M] by the addition of the 

damping matrix [D] given by Eq. 10c. The other necessary change is 

the addition of the Matrix [M - D (At/2)] as a product with [Aq*]^. . 

, i 

Thus since 
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'J-f 1 = 



[Mi] 



-1 



[m 2 ] 

0 



-1 G 



( 20 ) 



*w J 



Eq. 18 can be expressed in the form 



{Aq*} rt = [M + D r (At/2)] _1 f [M p - D r (At/2)] {Aq*} rt 
i+1 \ i 



( 21 ) 



+ ( C } 



rt.y ’ 



r = 1, 2,...R 
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C. Program Changes and Modification Logic 



The following routines of the IBM version of BR-1 have been modified 
for BR-1HR: MAIN, MEMBER, MTERM, QPLATE, ST0RE, DELTT, DEFLX, REVIV1 
and REVIV2. Two new subroutines were created: DPMASS and ADDAMP. The 

core size was increased from 250k bytes to 290k bytes. 

The flow of the logic of the modifications is as follows: 

1. Compute [M] in ST0RE (no change) 

2. Compute [D] in ST0RE 

3. Compute [M] 1 in MTERM (no change) 

4. Compute maximum time interval for numerical stability in DELTT 

based on [M] (no change) 

5. Take the inverse of [M]^ 1 to get [M]^ in DPMASS using INVS 

6. Compute [M + D y (At/2)] and [M - D r (At/2)] in DPMASS 

7. Compute [M r + D f (At/2)]" 1 in DPMASS using INVS 

8. Compute [M + D r (At/2)]" 1 [M - D r (At/2)] in DPMASS 

9. Compute [M +D (At/A)] " 1 [M -D (At/A)] {Aq*} in ADDAMP 

1C 1C 1C 1C ^ 

10. Compute j^Aq*}^ using Eq. 21 in DEFLX (no change) 

in 

The phrase "no change" means that the original procedure was used. 

When no damping is considered the modifications and additions are 
bypassed. 
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III. USER'S INSTRUCTIONS FOR BR-1HR 



The instructions for the use of the BR-1 code are given in Ref. 2. 
All of the instructions contained in that volume also apply to the modi- 
fied program BR-1HR. The time step for numerical stability of BR-1HR 
is identical to that of BR-1. The additional instructions required to 
use BR-1HR are as follows : 



1 . 



2 . 



Problem Control Card (page 4, Ref. 2) 

IHR (l5» Col. 66-70) - IHR = 0, no fluid is involved; the 

original BR-1 code is used. IHR = 1, 

(follows IREV) 

fluid is involved, the modifications 



are 

Rectangular Panel Card (page 8 
RH0CF (E8.4, Col. 55-62) - 

(between RH0 and Table 
C0DE) 



used. 

, Ref. 2) 

RH0CF is the product of y^., the 
fluid specific weight , and c, the 

sonic velocity of the fluid. The 

c -2 -1 

units of Y|i are lb^ -in. - sec. 

If the panel is not in contact 



with the fluid, RH0CF = 0. 



IV. SAMPLE PROBLEM 



A simply supported square plate is subjected to a step pressure load 
of the form 

p = o t = o 

t > o (22) 



ttX 



p = P sin — sin 

a a 



in 



Due to symmetry, only one quarter of the plate is considered: The 



parameters of the problem are: 
£ 

E = 10. 4 x 10 psi 

v = O.O 965 lb/in? 
u = 1/3 
h = 0.1 in. 
a =20 in. 



- Young's modulus 

- specific weight of the plate 

- Poisson’s ratio 

- thickness 

- length and width 



P = 0.01 lb/inf 

The load is sufficiently small such that the nonlinear effects are 
negligible. The plate is modeled with four elements as shown in Fig. 1. 

The equation governing the damped motion of the plate corresponding 
to Eq. 16 is 



Dy w + ^ — w + pew = Psin ^ sin ^ 
g y a a 



(23) 



where 



D = 



Eh~ 



12(l- v 2 ) 



V = 






4 4 

+ 2 rV 5 + 3* 

ax ay ay 



and g is the local acceleration due to gravity. The solution to Eq. 
23 is 
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SIMPLY SUPPORTED EDGE 



SIMPLY SUPPORTED EDGE 

2 3 



8 



I 



LINE OF SYMMETRY 
lOin. 



E - 1 0.4 X I O 6 psi u * l/ 3 

0.0965 Ib/in^ h=O.I In. 



FIGURE I. SAMPLE PROBLEM 
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LINE OF SYMMETRY 



W = w st [1 - e' C(ot cos (“\/l - C 2 ' (ut + cp) / cos <p} 



(24) 



where 



tan fp = - q/ 






a P ttX ny 

w , = r sin JJ — sin — 

st 4 a a 



4D 



tt 



r = 

2yh 




when the plate is initially at rest. The displacement at the center of 
the plate given by Eq. 24 is plotted in Figures 2, 3 } and 4 as a function 
of time for Q = 0, 0.666 and 270 corresponding to zero damping, less than 
critical damping and very heavy damping respectively. The corresponding 



2 

values of gpc for the fluid are 0, 4, and 1620 lb^ /(in -sec). Also 
plotted in Figs. 2-4 are the results from BR-1HR. The input data sheets 
and the print of the input data are given in Fig. 5- The execution time 
on the IBM 360/67, FORTRAN IV - Level H, was 8 min. and 26 sec. for 200 
time steps with Q = 270. The run with damping not considered took 
essentially the same length of time. 
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EXACT SOLUTION 




(c-OI X Nl) 



M 



FIGURE 2. TRANSVERSE DISPLACEMENT AT NODE 9 VERSUS TIME, f-0 (NO DAMPING) 



3. 




t (MSEC) 



FIGURE 3. TRANSVERSE DISPLACEMENT AT NODE 9 VERSUS 
TIME, 0.666 
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(g-OI X Nl ) 



1.4 




* (MSEC) 

FIGURE 4. TRANSVERSE DISPLACEMENT AT NODE 9 VERSUS 
TIME , J = 270 
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V. SUMMARY AND CONCLUSIONS 



The finite element digital computer code BR-1 developed by the Northrop 
Corporation for predicting the effects of internal air blast on typical 
combat aircraft skin- rib- stringer structures has been modified to include 
the effect of compressible fluid- structure interaction. The fluid- structure 
interaction is approximated by the piston theory wherein the effect of the 
fluid upon the structure is accounted for by introducing damping to the 
equations of motion of the structure. The modified code is called BR-1HR. 
This code, in conjunction with the NWC code for predicting hydraulic ram 
pressures, can be used to predict the structural response of aircraft fuel 
tanks subjected to penetrating bullets and fragments. 

All of the features of BR-1 exist in BR-1HR, and only two additional 
numbers are required for the input data. The modified code is operational 
on the IBM 360/67 in FORTRAN IV, level H, and requires 290K bytes of 
storage. A sample problem was executed to demonstrate the validity of the 
modified code for zero damping, less than critical damping, and very heavy 
damping. 
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